% Subroutine for computing stationary equilibrium wage
% need to set parameters by setparams.m
% October 2022
% Makoto Nirei

wh=1;
wl=0;
w=(wh+wl)/2; w_gap=1;
q=zeros(1,length(a)); tmp_w=zeros(1,length(a)); pb=q; pt=q; zt=q; qu=q;
while abs(w_gap)>tol_w
    for i=1:length(a)
        q(i)=fzero(@(x) findq(x,rm,delta,eta,pi,a(i),w),1.3);
        tmp_w(i)=vp(q(i),1-eta+mu/pi)/vp(q(i),mu/pi)...
            *(a(i)*vp(q(i),rm/pi+1-eta)/vp(q(i),rm/pi-eta))^(eta-1);
    end
    w_old=w;
    w_new = (eta-1)/eta * (mean(tmp_w))^(1/(eta-1));
    w_gap=w_new-w_old;
    if w_gap>0
        wl=w_old;
    else
        wh=w_old;
    end
    w=(wl+wh)/2;
end
for i=1:length(a)
    q(i)=fzero(@(x) findq(x,rm,delta,eta,pi,a(i),w),1.3);
    pb(i) = ((vp(q(i),rm/pi+1-eta)/vp(q(i),rm/pi-eta)*vp(q(i),-eta)...
        -vp(q(i),1-eta))*(eta-1)/rm/delta)^(1/(eta-1));
    pt(i) = q(i)*pb(i);
    zt(i) = pt(i)^(1-eta)-(w/a(i))*pt(i)^(-eta);
end
vt(1) = 0.5*((zt(1)+zt(2))/rho + (zt(1)-zt(2))/(rho+2*mu*zeta));
vt(2) = 0.5*((zt(1)+zt(2))/rho - (zt(1)-zt(2))/(rho+2*mu*zeta));

const=((1-eta)/(rm/pi+1-eta) *pt.^(rm/pi+1-eta) + eta/(rm/pi-eta)*(pt.^(-eta+rm/pi))./a *w)/rm;

% price cut
v_func = @(x,v,vj,c,eta,rm,pi,w,a,zeta,mu) ...
    c*x.^(-rm/pi)+x.^(1-eta)/(rm-pi*(eta-1))-x.^(-eta)*w/a/(rm-pi*eta)...
    +(zeta*vj+(1-zeta)*v)*mu/rm;
v_peak(1) = v_func(pt(1),vt(1),vt(2),const(1),eta,rm,pi,w,a(1),zeta,mu);
v_peak(2) = v_func(pt(2),vt(2),vt(1),const(2),eta,rm,pi,w,a(2),zeta,mu);
v_gap = @(x,v,vj,c,eta,rm,pi,w,a,zeta,mu, v_peak, delta_u)...
    v_func(x,v,vj,c,eta,rm,pi,w,a,zeta,mu)-v_peak+delta_u;
pu(1)=fzero(@(x) v_gap(x,vt(1),vt(2),const(1),eta,rm,pi,w,a(1),zeta,mu, v_peak(1), delta_u),[pt(1) pt(1)+1]);
pu(2)=fzero(@(x) v_gap(x,vt(2),vt(1),const(2),eta,rm,pi,w,a(2),zeta,mu, v_peak(2), delta_u),[pt(2) pt(2)+1]);
for i=1:length(a)
    qu(i) = pt(i)/pu(i);
end

% utility
ep_eta_a=zeros(1,length(a));
for i=1:length(a)
    ep_eta_a(i)=(pt(i)^(mu/pi-eta)-pb(i)^(mu/pi-eta))/(mu/pi - eta) / pb(i)^(mu/pi) /vp(q(i),mu/pi);
end
ss.w=w;
ss.ep=mean(ep_eta_a ./ a); %E[p^(-eta)/a]
ss.lambda= (mu/(q(1)^(mu/pi)-1) + mu/(q(2)^(mu/pi)-1))/2;
ss.MenuCost = ss.lambda *delta;
ss.CpN = (1- ss.MenuCost)/ ss.ep;
%alpha_U =1/(1+2*w/CpN); % calibrated value of alphaU when N=1/3. When pi=0.03, the calibragted value is 0.4213
ss.N = 1/(1+(1/alpha_U -1)* ss.CpN/ss.w);
ss.C = ss.N * ss.CpN;
% nominal return R is set to pi. So, rho-R+pi=rho.
ss.m = 1/(ss.C^(alpha_U*(1-gamma_U)-1)*(1-ss.N)^((1-alpha_U)*(1-gamma_U))*alpha_U*rho);
ss.q = q;
ss.pb=pb;
ss.pt=pt;
ss.pu=pu;